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Abstract. Observations of the cosmological large-scale structure provide well- 
established neutrino mass limits. We extend this argument to thermal relic axions. 
We calculate the axion thermal freeze-out temperature and thus their cosmological 
abundance on the basis of their interaction with pions. For hadronic axions we find 
a new mass limit ttIq < 1.05 eV (95% CL), corresponding to a limit on the axion 
decay constant of fa > 5.7 x 10^ GeV. For other models this constraint is significantly 
weakened only if the axion-pion coupling is strongly suppressed. For comparison we 
note that the same approach leads to J2 < 0-65 eV (95% CL) for neutrinos. 
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1. Introduction 

Low-mass cosmic relic particles provide a hot-dark matter component of the universe. 
After they decouple in the early universe, they freely stream and decrease the spectral 
power of primordial density fluctuations on small scales. This well-studied effect has 
been established as a standard method to derive neutrino mass limits where the results of 
different authors differ primarily in the used observational data sets and assumed priors 
on some of the cosmological parameters (for a recent review see Ref. [1]). Two of us have 
recently extended this method to several generic cases of other low-mass thermal relics 
that may have been produced in the early universe [2], for example weakly interacting 
bosons such as the hypothetical axions. 

For a generic low-mass particle species X that was once in thermal equilibrium, the 
crucial parameters that enter the structure-formation argument are the assumed particle 
mass mx, the number gx of internal degrees of freedom, the relevant quantum statistics 
(boson vs. fermion), and finally the number g^* of effective cosmic thermal degrees of 
freedom at the epoch when the X particles thermally decouple. These parameters 
imply the present-day cosmic mass density fix in these particles as well as their velocity 
distribution and primordial free-streaming scale. 

In Ref. [2] two of us considered thermal relic axions as a specific example and found 
an upper mass limit of 2-3 eV. The main purpose of the present paper is to sharpen 
this result in two ways. First, we use an updated set of observational data. Second, we 
take a closer look at the axion freeze-out in the early universe to obtain a more accurate 
value for in terms of the appropriate coupling constants. 

We begin in Sec. 2 with a brief summary of the relevant aspects of axion physics 
and continue in Sec. 3 with a determination of the freeze-out conditions as a function of 
the relevant axion coupling constants. In Sec. 4 we perform a likelihood analysis based 
on recent cosmological data to determine the allowed range of axion parameters. We 
finish in Sec. 5 with a summary and interpretation of our results. 

2. Axions 

Quantum chromodynamics is a CP-violating theory, implying that the neutron should 
have a large electric dipole moment, in conflict with the opposite experimental evidence. 
The most elegant solution of this "strong CP problem" was proposed by Peccei and 
Quinn (PQ) who showed that CP conservation is dynamically restored in the presence 
of a new global U(l) symmetry that is spontaneously broken at some large energy 
scale [3,4]. Weinberg [5] and Wilczek [6] realized that an inevitable consequence of the 
PQ mechanism is the existence of a new pseudoscalar boson, the axion, which is the 
Nambu-Goldstone boson of the PQ symmetry. This symmetry is explicitly broken at 
low energies by instanton effects so that the axion acquires a small mass. Unless there 
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are non-QCD contributions, perhaps from Planck-scale physics [7, 8] , the mass is 

_ z'/^ Um^ _ 6.0 eV 
1 + z fa /a/106 GeV 
where the energy scale is the axion decay constant or PQ scale that governs all axion 
properties, /tt = 93 MeV is the pion decay constant, and z — mu/m,,i is the mass ratio 
of the up and down quarks. We will follow the previous axion literature and assume a 
value z = 0.56 [9, 10], but we note that it could vary in the range 0.3-0.7 [11]. 

The PQ scale is constrained by various experiments and astrophysical arguments 
that involve processes where axions interact with photons, electrons, and hadrons. The 
interaction strength with these particles scales as f~^, apart from model-dependent 
numerical factors. If axions indeed exist, experimental and astrophysical limits suggest 
that /a ^0.6 X 10^ GeV and m„^0.01 eV [11, 12]. 

Perhaps the most robust limits are those on the axion interaction strength with 
photons and electrons because the properties of ordinary stars can be used to test 
anomalous energy losses caused by processes such as 7 + e~ — e~ + a or the Primakoff 
process 7 + Ze ^ Ze + a [11, 12]. These limits are so restrictive that processes involving 
the electron and photon couplings will not be significant for the thermalization of axions 
in the early universe in the post-QCD epoch. On the other hand, these interactions are 
also the most model dependent. Axion-models can be constructed where the axion- 
photon interaction is arbitrarily small. In "hadronic axion models" such as the KSVZ 
model [13, 14] there is no tree-level interaction with ordinary quarks and leptons, and 
in any case, the axion-electron coupling can be very small even in non-hadronic models 
such as the DFSZ model [15, 16]. Therefore, the restrictive stellar limits on the electron 
and photon couplings do not conclusively rule out axions with a PQ scale in the 10^ GeV 
range that is of interest to us. 

The axion-nucleon interaction strength is constrained by two different limits based 
on the observed neutrino signal of supernova 1987A [11,12]. The axionic energy loss 
caused by processes such as NN — > NNa excludes a window of the axion-nucleon 
interaction strength where the axionic contribution to the loss or transfer of energy 
would have been comparable to or larger than that of neutrinos. For a sufficiently 
large interaction strength, axions no longer compete with neutrinos for the overall 
energy transfer in the supernova core, but then would cause too many events in the 
water Cherenkov detectors that observed the neutrino signal. However, there is a 
narrow intermediate range of couplings, corresponding to fa around 10^ GeV, i.e. to an 
axion mass of a few eV, where neither argument is conclusive. In this "hadronic axion 
window," these elusive particles may still be allowed. This observation has previously 
led to the speculation that axions with eV-masses could contribute a significant hot 
dark-matter fraction to the universe [17]. 
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Assuming that axions do not couple to charged leptons, the main thermahzation 
processes in the post-QCD epoch involve hadrons [18] 

a + TT <-> TT + TT , (2) 

a + N ^ N + 7r . (3) 

Due to the scarcity of nucleons relative to pions, the pionic processes are actually by far 
the most important. Therefore, we need the axion-pion interaction which is given by a 
Lagrangian of the form [18] 

Can = + 7r°7r-a^7r+ - 27r+7r-a^7rO) d^a . (4) 

In J a 

In hadronic axion models where the ordinary quarks and leptons do not carry PQ 
charges, the couphng constant is [18] 

1 — z 

^"'^ = soTzj ■ 

In non-hadronic models an additional term enters that could, in principle, reduce Can- 
For the DFSZ model, the relevant a-7r'^-mixing strength was derived by Carena and 
Pcccei [19]. In the DFSZ model, however, axions in our /„ range are already excluded 
so that we mainly focus on hadronic models. 



3. Axion decoupling 
3.1. Pionic process 

Axions are an often-cited cold dark matter candidate if they interact so weakly (if /„ 
is so large) that they never achieve thermal equihbrium [23-27]. In that case coherent 
oscillations of the axion field are excited around the QCD transition when the PQ 
symmetry is explicitly broken by instanton effects. However, if fa is sufficiently small 
so that axions interact sufficiently strongly, they will thermalize and survive as thermal 
relics in analogy to neutrinos [18,20,21]. In the eV-range of axion masses that is of 
interest to us, thermahzation will happen after the QCD transition so that one has 
to consider axion interactions with hadrons rather than interactions with quarks and 
gluons that would be relevant at earlier epochs. 

Axions freeze out when their rate of interaction becomes slow compared to the 
cosmic expansion rate. As a criterion for axion thermal decoupling we use 

{Tah = H{T) , (6) 

where (^a)T is the axion absorption rate, averaged over a thermal distribution at 
temperature T, whereas H{T) is the Hubble expansion parameter at the cosmic 
temperature T. Our freeze-out criterion is accurate up to a constant of order unity. 
In principle, the full system of Boltzmann equations should be followed for all species, 
but the use of H{T) — (ra)^ only introduces an error of 10-20% in the decouphng 
temperature, consistent with the other approximations we will use. 
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As explained in Sec. 2, by far the most important process for axion thermalization 
is the pionic reaction of Eq. (2) which is of the form 1 + 2^3 + 4 with particle 1 the 
axion and particles 2-4 the pions. The average absorption rate is 

1 r (i^pi d^p2 d^Ps d^p4 
^ J {2nf2Ei (27r)32E2 {2tiY2E^ (27r)32^4 



X {2nf 5\p, + P2 " P3 - Pa) fMl + h) (1 + U) 



(7) 



where fj are the thermal occupation numbers. The axion number density in thermal 
equilibrium is na = / (i^pi/i/(27r)^ = (C3/7r^)T^. 

The interaction Eq. (4) implies the three pionic processes a + ti^ 7r+ + tt^, 
a+TT"*" — >• 7r"'" + 7r°, and a + 7r~ 7r~ + 7r°. Assuming equal mass ttIt^ for both charged and 
neutral pions, we find for the squared matrix element, summed over all three channels, 

E I-MI' = I i^] ' + t' + u'- Smi) , (8) 



where s = {pi + P2)^, t = {pi — Ps)^, and u = {pi — p^Y- 

We have performed the phase-space integration in Eq. (7) along the lines of the 
method described in Ref. [22]. Chang and Choi [18] give an explicit three-dimensional 
integral expression for the average rate of tttt — > 7ra when the final-state stimulation 
factors (l + /3)(l + /4) are neglected. Our numerical integration agrees with their result 
if we consider the same quantity, but in the following we use our full expression for the 
average absorption rate. Prom dimensional considerations one finds 

{Ta)T = A T'h{mjT), (9) 



, fafiT , 

where numerically A — 0.215. The function is normalized to /i(0) = 1 (Fig. 1). 




9 10 



Figure 1. Dimensionless axion absorption rate /i(/x) as defined by Eq. (9). 
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3.2. Expansion rate 



The cosmic expansion rate is given by the Priedmann equation as — SttGnp/S with 
Gn the Newton constant and p the mass-energy density. For the radiation dominated 
epoch one usually writes [28] 

d'p (m2 + p2)V2 



TT 



E 9j 

j=particles 



/ 



exp 



(m2 + p2)V2/r 



± 1 



(10) 



where g*{T) denotes the effective number of thermal degrees of freedom that are excited 
at the epoch with temperature T. Further, gj is the number of internal degrees of 
freedom of a given particle species j with mass mj. The ±1 in the denominator applies 
to bosons and fermions, respectively. We assume vanishing chemical potentials for all 
particles. For the conditions of interest, the particle-antiparticle asymmetry is small 
even for nucleons. The Friedmann equation is thus 

1/2 



H 



47^ 

45 



9*{T) 



rp-2 

mpi 



where mpi = G^^'^ is the Planck mass. 

We consider all particles with masses up to the nucleon mass as hsted in Table 1. In 
the upper panel of Fig. 2 we show the function g^(T), assuming that only the particles of 
Table 1 contribute and ignoring the color deconfinement transition that probably takes 
place at a temperature somewhat below 200 MeV. 



Table 1. Particles that contribute to 5* during the post-QCD epoch 



Particles 


Mass [MeV] 


Spin 


9j 


7 





1 


2 







1 
2 


6 


e± 


0.511 


1 

2 


4 




106 


1 

2 


4 


ttO, ^± 


138 





3 


i^o, K\ 


494 





4 


V 


547 





1 




771 


1 


9 


CO 


782 


1 


3 


K*°, K*°, K*^ 


892 


1 


12 


n, p. n. J) 


9.'-i8 


J_ 


8 



The effective number of thermal degrees of freedom g^ that appears in the 
Friedmann equation is the quantity that governs the expansion rate. However, if axions 
freeze out at a certain epoch, their number density at late times is governed by g^s, the 
effective number of thermal degrees characterizing the entropy at the freeze-out epoch. 
The entropy density of an ideal gas is given by the general expression {p+P)/T oc g^sT^ 
with P the pressure. The pressure is given by the same integral expression as the energy 
density in Eq. (10) after replacing Ej = (m^ -|- p^)^/^ by | p • v = ^p'^/Ej. Even though 
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Figure 2. Effective number of thermal degrees of freedom in the early universe during 
the post-QCD epoch, assuming the particle content listed in Table 1. Upper panel: g^. 
Lower Panel: g*s/g*- 



we assume that all particles are in thermal equilibrium at the same temperature, there 
will be a difference between g^, and (7*5 because some of the contributing particles are not 
massless. In the lower panel of Fig. 2 we show the ratio g^^s/ 9*- Since the deviation of g^, 
from g^s is at most a few percent for the conditions of interest, we will henceforth ignore 
the difference between the two quantities and always use g^. Moreover, since axions 
themselves contribute only a single degree of freedom we neglect their contribution to 
g^ for simplicity. 



3.3. Freeze-out conditions 

We now combine our result for the cosmic expansion rate in the post-QCD epoch 
with that for the axion absorption rate and determine the freeze-out conditions from 
Eq. (6). As an example we show H{T) and (ya)T in Fig. 3, assuming a PQ scale of 
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fa — 10^ GeV and the hadronic axion-pion coupling of Eq. (5). From the intersection 
point we determine Td, the axion decouphng temperature, and the corresponding g^. 
We show these quantities as functions of fa in Fig. 4 and in Table 2. 

Moreover, from g^, at decoupling one can determine the present-day number density 
of axions by virtue of the relation 

5(*5 (today) 



n„, = 



5f*s(decoupling) 



(12) 



where = 411 cm~^ [11] is the present-day density of cosmic microwave photons 
and today g^^^s — 3.91 [28]. We show Ua as a function of in Fig. 4 and give some 
numerical values in Table 2. For comparison, we note that the present-day neutrino 
number density, determined by standard neutrino freeze-out, is rii, ~ 112 cm~^. 




20 40 



60 80 100 120 140 160 180 200 
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Figure 3. Average axion absorption rate for /„ = 10'' GeV and cosmic expansion rate 
as a function of the cosmic temperature. Rates are in units of T^/mpi. 



Table 2. Thermal decouphng conditions for axions, assuming the hadronic axion-pion 
couphng of Eq. (5). Ua is the expected present-day axion density. 



fa [GeV] Td [MeV] 5,(Td) Ua [cm-^j 



3 X 10^ 


13.37 


10.84 


74.14 


1 X 10^ 


15.30 


10.93 


73.50 


3 X 10" 


17.63 


11.10 


72.39 


1 X 10^ 


21.21 


11.46 


70.11 


3 X 10^ 


26.06 


12.06 


66.63 


1 X 10^ 


34.75 


13.15 


61.08 


3 X 10*^ 


49.12 


14.54 


55.24 


1 X 10^ 


81.61 


16.43 


48.88 


3 X 10^ 


145.31 


21.10 


38.08 
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Figure 4. Axion freeze-out, assuming the a-ir coupling of Eq. (5). Top: Decoupling 
temperature. Middle: Corresponding (7* . Bottom: Present-day axion density. 
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4. Likelihood analysis 

4.1. Theoretical predictions 

In order to derive limits on the axion parameters we compare theoretical power spectra 
for the matter distribution and the CMB temperature fluctuations with observational 
data in analogy to previous works by one of us [29-31]. The predicted power spectra 
arc calculated with the publicly available CMBFAST package [32]. As in our previous 
paper on general low-mass thermal relics [2] we have modified the code to include a 
scalar boson that decouples earlier than neutrinos, i.e. we assume a number density and 
velocity distribution according to the freeze-out epoch. 

As a set of cosmological parameters apart from the axion characteristics we 
choose the matter density Qu, the baryon density Qb, the Hubble parameter Ho, 
the scalar spectral index of the primordial fluctuation spectrum Ug, the optical depth 
to rcionization r, the normalization of the CMB power spectrum Q, and the bias 
parameter b. We restrict our analysis to geometrically flat models Q = Qu + = 1- 
The cosmological parameters and their assumed priors are listed in Table 3. The 
actual marginalization over parameters was performed using a simulated annealing 
procedure [33]. The cosmological parameters in our model corresponds to the simplest 
ACDM model which fits present data. While our results would be changed slightly by 
including additional parameters, there would be no significant changes [2]. 

Likelihoods are calculated from so that for 1 parameter estimates, 68% confidence 
regions are determined by Ax^ — ~ Xo — ^5% regions by Ax^ — 4. Here, Xo 

refers to the best-fit model. 

Table 3. Priors on cosmological parameters used in the likelihood analysis. 



Parameter Prior Distribution 





FOa 1 


Fixed 




0.1-1 


Top Hat 


h 


0.72 ±0.08 


Gaussian 




0.014-0.040 


Top hat 


Us 


0.6-1.4 


Top hat 


T 


0-1 


Top hat 


Q 




Free 


b 




Free 



4-2. Cosmological Data 

4-2.1. Large Scale Structure. At present there are two large galaxy surveys of 
comparable size, the Sloan Digital Sky Survey (SDSS) [34, 35] and the 2 degree Field 
Galaxy Redshift Survey (2dFGRS) [36]. Once the SDSS is completed it will be 
significantly larger and more accurate than the 2dFGRS. We will only use data from 
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SDSS, but the results would be almost identical with 2dFGRS data. We use only data 
on scales larger than k — 0.15h Mpc~^ to avoid problems with non-linearity. 

4-2.2. Cosmic Microwave Background. The CMB temperature fluctuations are 
conveniently described in terms of the spherical harmonics power spectrum Cf^ = 
{\0'im\'^), where ^{9,(j)) = Y,im (^imyimif^ , 4>) ■ Since Thomson scattering polarizes light, 
there are also power spectra coming from the polarization. The polarization can be 
divided into a curl-free (E) and a curl (B) component, yielding four independent power 
spectra: Cf^, Cf ^, Cf ^, and the T-E cross-correlation Cf^. 

The WMAP experiment has reported data only on Cf^^ and Cj^^ as described in 
Refs. [37-41]. We have performed our likelihood analysis using the prescription given by 
the WMAP collaboration [38-41] which includes the correlation between different C;'s. 
Foreground contamination has already been subtracted from their published data. 

4.2.3. Supernova luminosity distances. We use the "gold" data set compiled and 
described by Riess et al. [42] consisting of 157 SNe la using a modified version of the 
SNOG package [43]. 

4 .2. 4- The Lyman-a forest. We have furthermore used the Lyman-cc forest data from 
Croft et al. [44]. The error bars on the last three data points have been increased in 
the same fashion as was done by the WMAP collaboration [37-41], in order to make 
them compatible with the analysis of Gnedin and Hamilton [45]. While there is some 
controversy about the use of Lyman-o; data for parameter estimation it should be noted 
that their inclusion has only a very small quantitative effect on our results as will become 
clear below. 

4-2.5. Hubble parameter. Additionally we use data on the Hubble parameter from the 
HST key project [46], h = //o/(100 km s"^ Mpc"^) = 0.72 ± 0.08. 

4-3- Limits on axion parameters 

After marginahzing over the cosmological parameters shown in Table 3, the 68% and 
95% CL allowed regions of axion parameters are shown in Fig. 5. In the upper panel 
we have included the full data set described above, while in the lower panel we have 
removed the Lyman-a data that is perhaps our most uncertain input. 

The Lyman-a data only has a significant impact for high values of because 
it measures very small scales. The large-scale structure data loses sensitivity at 
k^O.lbh Mpc~^, whereas the Lyman-a data probes scales which are about an order 
of magnitude smaller. For high g^,, the particle mass for a given value of ^ah^ is higher 
and therefore the free-streaming length is correspondingly smaller. The change in power 
spectrum amplitude occurs at A; > 0.15 /i Mpc~^ and is only measurable in the Lyman-a 
data when g^ ^ 60-70. 




0.004 0.006 0.008 0.010 0.012 0.014 




0.004 0.006 0.008 0.010 0.012 0.014 

Figure 5. Likelihood contours for the allowed axion parameters. Everything to the 
right of the dark shaded region is excluded at the 68% CL, and everything to the right 
of the light shaded region is excluded at the 95% CL. Contours of equal axion mass 
with the indicated TOq values in eV are shown as dashed lines. Upper panel: All data 
included. Lower panel: Lyman-a data excluded. 

In Fig. 6 we show the same analysis as in the upper panel of Fig. 5, but now 
transformed to the ma-(7*-plane. The hadronic axion model is shown as a thick solid 
line. It lies at relatively small values of for the relevant mass range, meaning that 
there will be very little quantitative difference if the Lyman-a data is not used. 
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12 3 4 

m„ (eV) 

Figure 6. Likelihood contours including all data, equivalent to the upper panel of 
Fig. 5, transformed to the parameter plane of and m^. The relation between and 
nia for hadronic axions is shown as a thick solid line. 




0.5 1.0 1.5 2.0 2.5 3.0 
m„ (eV) 

Figure 7. Likelihood contours including all data, equivalent to the upper panel of 
Fig. 5, transformed to the parameter plane of fa and ma- We have used the connection 
between g^, and fa for hadronic axions as shown in Fig. 4 and Table 2. The standard 
Wo-/a-relation of Eq. (1) is shown as a thick solid hue. 
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Likewise, in Fig. 7 we show the analysis in the m(j-/(j-plane, where the relationship 
between and g^, shown in Fig. 4 and Table 2 has been used. In this figure, the true 
independent variable on the vertical axis is the axion-pion coupling that was transformed 
to a value for the PQ scale assuming the hadronic case Eq. (5). For other models with 
different Ca-K-, the vertical axis in Fig. 7 must be re-scaled accordingly. Moreover, the 
standard relationship between ma and fa of Eq. (1) and the hadronic axion-pion coupling 
of Eq. (5) put axions on the thick solid line. 

In this case, where the most specific assumptions about the underlying axion model 
were made, the axion parameter space collapses to one dimension, leading to a more 
restrictive mass limit. In Fig. 8 we show a one-dimensional likelihood analysis for this 
case, i.e. for the hadronic axion model. As a result one finds a 95% CL excluded region 
for nia of 

ma < 1.05 eV (95% CL), (13) 

based on all the available observational data. For comparison we note that the same 
data and the same analysis method provide 

^m^<0.65eV (95% CL) (14) 

for neutrinos [31]. 

12 
10 
8 

01 

X 

< 6 

4 

2 


0.0 0.5 1.0 1.5 2.0 2.5 
m„ (eV) 

Figure 8. Onc-dimcnsional likelihood analysis for the hadronic axion model obeying 
Eq. (1). Shown is Ax^ = — as a function of where Xo refers to the best-fit 
model. In this ID analysis the 95% excluded region is at Ax^ > 4, and the 68% 
excluded region at Ax^ > 1, as indicated by the horizontal dashed lines. 
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5. Summary 

We have studied structure-formation limits on axions or other scalar particles that have 
a small mass and that couple to pions by virtue of the Lagrangian Eq. (4). Our main 
result is the exclusion range in the two-dimensional fa-ma parameter space shown in 
Fig. 7 where the axion-pion coupling strength of Eq. (5) was assumed, i.e. this is an 
exclusion range in the parameter space of axion mass and axion-pion coupling. Assuming 
in addition the standard rria-fa relationship of Eq. (1), wc find a limit on the axion mass 
of m„ < 1.05 eV (95% CL), corresponding to /„ > 5.7 x 10^ GeV. 

A comparable axion mass limit is obtained from the requirement that excessive 
energy losses of horizontal branch stars in globular clusters should be avoided [12]. 
However, the globular-cluster limit depends on the axion-photon interaction that is 
rather model dependent. 

The neutrino observations of supernova 1987A provide restrictive limits on the 
axion-nucleon interaction, suggesting ma ^ 0.01 eV or /„ ^ 0.6 x 10^ GeV if one assumes 
generic coupling strengths between axions and nucleons. However, the supernova 
argument may leave open the "hadronic axion window" of about 3 x 10^ GeV ^ /a ^ 3 x 
10^ GeV [17, 18]. Our result closes this window, at least for generic values of the axion- 
nucleon and axion-pion couplings. 

More importantly wc have provided a new limit on axion parameters based on 
a different interaction channel than previous limits and based on different data and 
assumptions. Every experimental measurement and every astrophysical or cosmological 
argument has its own systematic uncertainties and its own recognized or un-recognized 
loop holes. Therefore, to corner axions it is certainly important to use as many 
independent interaction channels and as many different approaches as possible. 

An ongoing experimental axion search, the GERN Axion Solar Telescope (CAST), is 
based on the axion-photon interaction. It has reported first limits for ^0.02 cV [47]. 
In the second phase it will cover axion masses up to approximately 1 cV. It is noteworthy 
that our new limit is almost identical with the upper mass range that can be reached 
with CAST, i.e. the CAST search and our new limit are nicely complementary. 
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